/*********************************************************************
*                                                                    *
* inidat -- provides initial data                                    *
*                                                                    *
* Parameters:                                                        *
* y_a_ptr        -- pointer where to write the initial data          *
* initial_time   -- initial time                                     *
*                                                                    *
* Returns: pointer to field_array where data was writen              *
*                                                                    *
*********************************************************************/

#include "first_macro_2d.h"  /* Where global parameters are defined */
#include "structs_2d.h"      /* Where structures are defined */
#include "derivs_2d.h"       /* Where derivatives functions are defined */
#include "gen_2d.h"
#ifdef SV
#include "java_ser.h"
#include "java_fser.h"
#endif



/* struct field_array *inidat(struct field_array *y_a_ptr) { */

void inidat(struct field_array *y_a_ptr, 
	    struct GRID_PAR *grid_ptr,
	    struct INI_PAR *ini_par_ptr){

/* -------> grid parameters <--------------------------*/


  
    int ni_1 = (*grid_ptr).start_grid_1; 
    int nf_1 = (*grid_ptr).final_grid_1;  
    int ni_2 = (*grid_ptr).start_grid_2;  
    int nf_2 = (*grid_ptr).final_grid_2; 

    int n_gridpts_1 = nf_1 - ni_1;
    int n_gridpts_2 = nf_2 - ni_2;

/*     int grid = (*grid_ptr).grid; */

#ifdef EXCISION
    int exc_n1 = (*grid_ptr).exc_grid_1;
    int exc_n2 = (*grid_ptr).exc_grid_2;
    int exc_n1i = (nf_1-ni_1-1)/2-exc_n1+1; /* take out a -1 for visual compatibility with the rest of code */
    int exc_n2i = (nf_2-ni_2-1)/2-exc_n2+1;
    int exc_n1f = (nf_1-ni_1-1)/2+exc_n1; 
    int exc_n2f = (nf_2-ni_2-1)/2+exc_n2;
#endif

    FLOAT x_i = (*grid_ptr).initial_x;
    FLOAT x_f = (*grid_ptr).final_x;
    FLOAT y_i = (*grid_ptr).initial_y;
    FLOAT y_f = (*grid_ptr).final_y;
    FLOAT swap, sigma_b02, r2;
    FLOAT g_up_xx, g_up_yy, g_up_xy, g_up_yx;
 
    int coorx = 1;
    int coory = 2;

/* -------> initial data parameters <------------------*/

 FLOAT twoPIdN1 = 2.*PI*(x_f-x_i)/(FLOAT)(nf_1-ni_1-1);
 FLOAT one_dN1 = (x_f-x_i)/(FLOAT)(nf_1-ni_1-1);


 FLOAT twoPIdN2 = 2.*PI*(y_f-y_i)/(FLOAT)(nf_2-ni_2-1);
 FLOAT one_dN2 = (y_f-y_i)/(FLOAT)(nf_2-ni_2-1);

 
  /* Parameters from main */


  FLOAT a0 = (*ini_par_ptr).a0;                
  FLOAT k_a_10 = (*ini_par_ptr).k_a_10; 
  FLOAT k_a_20 = (*ini_par_ptr).k_a_20;   
  FLOAT shift_a0 = (*ini_par_ptr).shift_a0;          

/* Amplitude sin(k_a1* x + shift_a1) in U1 */

  FLOAT a1 = (*ini_par_ptr).a1;                
  FLOAT k_a_11 = (*ini_par_ptr).k_a_11;  
  FLOAT k_a_21 = (*ini_par_ptr).k_a_21; 
  FLOAT shift_a1 = (*ini_par_ptr).shift_a1;          

/* Amplitude of cos(k_a0* x + shift_a0) in U0 */

  FLOAT c0 = (*ini_par_ptr).c0;                
  FLOAT k_c_10 = (*ini_par_ptr).k_c_10;    
  FLOAT k_c_20 = (*ini_par_ptr).k_c_20; 
  FLOAT shift_c0 = (*ini_par_ptr).shift_c0;          

/* Amplitude of cos(k_c1* x + shift_c0) in U1 */
  FLOAT c1 = (*ini_par_ptr).c1;                
  FLOAT k_c_11 = (*ini_par_ptr).k_c_11; 
  FLOAT k_c_21 = (*ini_par_ptr).k_c_21;
  FLOAT shift_c1 = (*ini_par_ptr).shift_c1;          

/* Amplitude of exp(cos(k_b0*x)^2/sigma_b0) in U0 */
  FLOAT b0 = (*ini_par_ptr).b0;                
  FLOAT sigma_b0 = (*ini_par_ptr).sigma_b0;
  FLOAT c0_1 = (*ini_par_ptr).c0_1;
  FLOAT c0_2 = (*ini_par_ptr).c0_2;

/* Amplitude of exp(cos(k_b1*x)^2/sigma_b1) in U1 */
  FLOAT b1 = (*ini_par_ptr).b1;                
  FLOAT sigma_b1 = (*ini_par_ptr).sigma_b1;
  FLOAT c1_1 = (*ini_par_ptr).c1_1;
  FLOAT c1_2 = (*ini_par_ptr).c1_2;

/* Global wave motion */

  FLOAT v1 = (*ini_par_ptr).v1;
  FLOAT v2 = (*ini_par_ptr).v2;

  int initial_data_type = (*ini_par_ptr).initial_data_type;


  /* --------------------------------------------------------------------------------------- */

#ifdef DEBUG_INIDAT
  printf("Begining of function inidat\n");
#endif

#ifdef EXCISION
{
if ((FLOAT)exc_n1*one_dN1 > 0.38*m) 
printf("excision region too large, R_exc_1 = %f \n", (FLOAT)exc_n1*one_dN1);

if ((FLOAT)exc_n2*one_dN2 > 0.38*m)
printf("excision region too large, R_exc_2 = %f \n", (FLOAT)exc_n2*one_dN2);

}
#endif

  /* --------------------------------------------------------------------------------------- */

  /* The cross-like coordinates */

switch ((*grid_ptr).grid)
{
 case 3:
  swap = x_i - (x_f-x_i);
  x_f = x_f - (x_f-x_i);
  x_i = swap;
  break;
 case 4:
  swap = y_i + (y_f-y_i);
  y_f = y_f + (y_f-y_i);
  y_i = swap;
  break;
 case 5:
  swap = y_i - (y_f-y_i);
  y_f = y_f - (y_f-y_i);
  y_i = swap;
  break;
 default:
  swap = x_i + (x_f-x_i)*(*grid_ptr).grid;
  x_f = x_f + (x_f-x_i)*(*grid_ptr).grid;
  x_i = swap;
  break;
}


/*---------> values for different fields <-------------*/


 

/* struct field_array y; */

/* first the time */

/* y.a.time = initial_time; */

(*y_a_ptr).time = (*grid_ptr).initial_time;



  
  /*----------- Initial data for the scalar field ----------- */

  switch (initial_data_type){
  case 0: // give data using cross-like coordinates with fourth order polynomial



      {register int g_ind1, g_ind2;
   FLOAT x=0.;
   FLOAT y=0.;
   FLOAT r=0.;

  for (g_ind1 = ni_1; g_ind1 < nf_1; ++g_ind1) {
    for (g_ind2 = ni_2; g_ind2 < nf_2; ++g_ind2) {
	  x =  x_i+(FLOAT)(g_ind1)*one_dN1;
	  y =  y_i+(FLOAT)(g_ind2)*one_dN2;

/*     (*y_a_ptr).u[M_A][IND_2(g_ind1,g_ind2)] = 	 */
/*        (a0*sin(  (FLOAT)g_ind1*twoPIdN1*k_a_10  */
/* 		 + (FLOAT)g_ind2*twoPIdN2*k_a_20  */
/* 		 + shift_a0) */
/*      + c0*cos(  (FLOAT)g_ind1*twoPIdN1*k_c_10  */
/* 		+ (FLOAT)g_ind2*twoPIdN2*k_c_20 */
/* 		+ shift_c0) */
/*        )*b0*exp(-sigma_b0*sigma_b0*( */
/* 		  (x-c0_1)*(x-c0_1) */
/* 		  + (y-c0_2)*(y-c0_2) */
/* 		  ) */
/* 		); */

	  r2 = (x-c0_1)*(x-c0_1) + (y-c0_2)*(y-c0_2);
	  sigma_b02 = sigma_b0*sigma_b0;
	 if (r2 <= sigma_b02){
	     (*y_a_ptr).u[M_A][IND_2(g_ind1,g_ind2)] = 
		 b0*pow((r2-sigma_b02),4)/pow(sigma_b02,4); 
	     (*y_a_ptr).u[M_B][IND_2(g_ind1,g_ind2)] = 
		 b0*pow((r2-sigma_b02),4)/pow(sigma_b02,4); 
	     (*y_a_ptr).u[M_C][IND_2(g_ind1,g_ind2)] = 
		 b0*pow((r2-sigma_b02),4)/pow(sigma_b02,4); 

	 }
	 else
	     {
		(*y_a_ptr).u[M_A][IND_2(g_ind1,g_ind2)] = 0.0;
		(*y_a_ptr).u[M_B][IND_2(g_ind1,g_ind2)] = 0.0;
		(*y_a_ptr).u[M_C][IND_2(g_ind1,g_ind2)] = 0.0;
	}
    }
  }
      }

    

    break;

  case 1: // cos(2phi) dependence 
      {
      

       FLOAT fz(FLOAT z);

      switch ((*grid_ptr).grid){ // patches 4 and 5 are the ones that have the poles
      case 0:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
        FLOAT z=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0+(FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0+(FLOAT)(grid_ind2)*one_dN2; 
              z = y/sqrt(1.+ x*x + y*y);
	      
				    
	      (*y_a_ptr).u[M_A]GRD_IND = (1.-x*x)/(1.+x*x)*cos(fz(z));
	      (*y_a_ptr).u[M_B]GRD_IND = 2.*x/(1.+x*x)*cos(fz(z));
	      (*y_a_ptr).u[M_C]GRD_IND = sin(fz(z));

	    }}}
	break;
				
      case 2:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
        FLOAT z=0.;
		
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2; 
              z = y/sqrt(1.+ x*x + y*y);
	      
				    
	      (*y_a_ptr).u[M_A]GRD_IND = (1.-x*x)/(1.+x*x)*cos(fz(z));
	      (*y_a_ptr).u[M_B]GRD_IND = 2.*x/(1.+x*x)*cos(fz(z));
	      (*y_a_ptr).u[M_C]GRD_IND = sin(fz(z));

	    }}}
	break;
		
      case 1:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
        FLOAT z=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =   -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =   -1.0 + (FLOAT)(grid_ind2)*one_dN2;      
              z = y/sqrt(1.+ x*x + y*y);
	      
				    
	      (*y_a_ptr).u[M_A]GRD_IND = -(1.-x*x)/(1.+x*x)*cos(fz(z));
	      (*y_a_ptr).u[M_B]GRD_IND = -2.*x/(1.+x*x)*cos(fz(z));
	      (*y_a_ptr).u[M_C]GRD_IND = sin(fz(z));

	    }}}
	break;
				
      case 3:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
        FLOAT z=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =   -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =   -1.0 + (FLOAT)(grid_ind2)*one_dN2;  
              z = y/sqrt(1.+ x*x + y*y);
	      
				    
	      (*y_a_ptr).u[M_A]GRD_IND = -(1.-x*x)/(1.+x*x)*cos(fz(z));
	      (*y_a_ptr).u[M_B]GRD_IND = -2.*x/(1.+x*x)*cos(fz(z));
	      (*y_a_ptr).u[M_C]GRD_IND = sin(fz(z));
			
	    }}}
	break;
		
      case 4:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
        FLOAT z=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =   -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =   -1.0 + (FLOAT)(grid_ind2)*one_dN2; 
              z = 1./sqrt(1.+x*x+y*y);
	     
				    
	if(x*x+y*y >= 0.00001){
	      (*y_a_ptr).u[M_A]GRD_IND = -(x*x-y*y)/(x*x+y*y)*cos(fz(z));
	      (*y_a_ptr).u[M_B]GRD_IND = -2.*x*y/(x*x+y*y)*cos(fz(z));
	      (*y_a_ptr).u[M_C]GRD_IND = sin(fz(z));
	}
	else{
	      (*y_a_ptr).u[M_A]GRD_IND = 0.0;
	      (*y_a_ptr).u[M_B]GRD_IND = 0.0;
	      (*y_a_ptr).u[M_C]GRD_IND = sin(fz(z));
	}														
	    }}}

	break;
      case 5:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
        FLOAT z=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =   -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =   -1.0 + (FLOAT)(grid_ind2)*one_dN2; 
              z = -1./sqrt(1.+x*x+y*y);
               
	
	if(x*x+y*y >= 0.00001){	
	      (*y_a_ptr).u[M_A]GRD_IND = -(x*x-y*y)/(x*x+y*y)*cos(fz(z));
	      (*y_a_ptr).u[M_B]GRD_IND = 2.*x*y/(x*x+y*y)*cos(fz(z));
	      (*y_a_ptr).u[M_C]GRD_IND = sin(fz(z));
	}
	else{
	      (*y_a_ptr).u[M_A]GRD_IND = 0.0;
	      (*y_a_ptr).u[M_B]GRD_IND = 0.0;
	      (*y_a_ptr).u[M_C]GRD_IND = sin(fz(z));
	}										
	    }}}
	break; 
      }
	}
    
    break;
    
  case 2: // give data using cross-like coordinates with 8th order polynomial 
      {register int g_ind1, g_ind2;
   FLOAT x=0.;
   FLOAT y=0.;
   FLOAT r=0.;

  for (g_ind1 = ni_1; g_ind1 < nf_1; ++g_ind1) {
    for (g_ind2 = ni_2; g_ind2 < nf_2; ++g_ind2) {
	  x =  x_i+(FLOAT)(g_ind1)*one_dN1;
	  y =  y_i+(FLOAT)(g_ind2)*one_dN2;

	  r2 = (x-c0_1)*(x-c0_1) + (y-c0_2)*(y-c0_2);
	  sigma_b02 = sigma_b0*sigma_b0;
	 if (r2 <= sigma_b02){
	     (*y_a_ptr).u[M_A][IND_2(g_ind1,g_ind2)] = 
		 b0*pow((r2-sigma_b02),8)/pow(sigma_b02,8); 
	 }
	 else
	     {(*y_a_ptr).u[M_A][IND_2(g_ind1,g_ind2)] = 0.0;}
    }
  }
      }
    break;

/* -------------------------------------------------------------*/
/* -------------------------------------------------------------*/
/* CALSE l=2,m-1 multipole */
/* -------------------------------------------------------------*/
/* -------------------------------------------------------------*/
  case 3: // l=2,m=1 multipole
      switch ((*grid_ptr).grid){ // patches 4 and 5 are the ones that have the poles
      case 0:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2; 
				    
	      (*y_a_ptr).u[M_A]GRD_IND = c0*(y/(1.0+x*x+y*y));
	    }}}
	break;
				
      case 2:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2; 
			
	      (*y_a_ptr).u[M_A]GRD_IND = c0*(-y/(1.0+x*x+y*y));
	    }}}
	break;
		
      case 1:
	{register int grid_ind1, grid_ind2, grid_ind3;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
				    
	      (*y_a_ptr).u[M_A]GRD_IND = c0*(-y*x/(1.0+x*x+y*y));
	    }}}
	break;
				
      case 3:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
			
	      (*y_a_ptr).u[M_A]GRD_IND = c0*((y*x)/(1.0+x*x+y*y));
	    }}}
	break;
		
      case 4:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
				    
	      (*y_a_ptr).u[M_A]GRD_IND = c0*(-y/(1.0+x*x+y*y));
	    }}}

	break;
      case 5:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
	
	      (*y_a_ptr).u[M_A]GRD_IND = c0*(-y/(1.0+x*x+y*y));
	    }}}
	break; 
      }
    
    break;


/* -------------------------------------------------------------*/
/* -------------------------------------------------------------*/

/* CALSE l=4 multipole */

  case 4: // l=4 multipole
    {register int grid_ind1, grid_ind2;
     FLOAT x, y, theta, phi;
//    printf("at inidat_type = 4, k_a_10 = %f, k_a_20 =%f \n",k_a_10,k_a_20);
    for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
      for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {

	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
	      theta = theta_cub(x,y,(*grid_ptr).grid) ;
          phi = phi_cub(x,y,(*grid_ptr).grid);

	  (*y_a_ptr).u[M_A]GRD_IND =  a0*sin(k_a_10*theta)*sin(k_a_20*phi + shift_a0);
	  (*y_a_ptr).u[M_B]GRD_IND =  a1*sin(k_a_11*theta)*cos(k_a_21*phi + shift_a1);
	  (*y_a_ptr).u[M_C]GRD_IND =  c0*cos(k_c_10*theta)*cos(k_c_20*phi + shift_c0);
//    printf("x=%f,y=%f,grid=%d,theta=%f,phi=%f\n",x,y,(*grid_ptr).grid,theta,phi);

	}}}
    break;

  default: //unknown initial data type
    printf("scalar field, unknown data type\n %d", initial_data_type);
    exit(1);
    break;
  }


// The space derivatives 




  /*--------- data for the x derivative of u[0] --------------*/


DERIV(grid_ptr,(struct field *)&(*y_a_ptr).u[M_A],(struct field *)&(*y_a_ptr).u[M_A_X],&coorx);

DERIV(grid_ptr,(struct field *)&(*y_a_ptr).u[M_B],(struct field *)&(*y_a_ptr).u[M_B_X],&coorx);

DERIV(grid_ptr,(struct field *)&(*y_a_ptr).u[M_C],(struct field *)&(*y_a_ptr).u[M_C_X],&coorx);

  /*---------- data for the y derivative of u[0] ------------- */
DERIV(grid_ptr,(struct field *)&(*y_a_ptr).u[M_A],(struct field *)&(*y_a_ptr).u[M_A_Y],&coory);

DERIV(grid_ptr,(struct field *)&(*y_a_ptr).u[M_B],(struct field *)&(*y_a_ptr).u[M_B_Y],&coory);

DERIV(grid_ptr,(struct field *)&(*y_a_ptr).u[M_C],(struct field *)&(*y_a_ptr).u[M_C_Y],&coory);





  /*----------- data for the time derivative of u[0] ---------------- */
  switch (initial_data_type){/* initial_data_type switch */

    case 0: {// give data using cross-like coordinates
      {register int grid_ind1, grid_ind2;
      FLOAT x=0.;
      FLOAT y=0.;

      for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {

	    x =  x_i+(FLOAT)(grid_ind1)*one_dN1;
	    y =  y_i+(FLOAT)(grid_ind2)*one_dN2; 

	  }}
      }}
      break;
  
    case 1: {// l=2 multipole
      switch ((*grid_ptr).grid){ // patches 4 and 5 are the ones that have the poles
      case 4:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
				    
	      	(*y_a_ptr).u[M_A_T]GRD_IND = c1*
			(-0.5 + 1.5/(1.0 + x*x + y*y));
		(*y_a_ptr).u[M_B_T]GRD_IND = c1*
			(-0.5 + 1.5/(1.0 + x*x + y*y));
		(*y_a_ptr).u[M_C_T]GRD_IND = c1*
			(-0.5 + 1.5/(1.0 + x*x + y*y));
	    }}}
	break;
				
      case 5:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
				    		    
	      	(*y_a_ptr).u[M_A_T]GRD_IND = c1*
			(-0.5 + 1.5/(1.0 + x*x + y*y));
		(*y_a_ptr).u[M_B_T]GRD_IND = c1*
			(-0.5 + 1.5/(1.0 + x*x + y*y));
		(*y_a_ptr).u[M_C_T]GRD_IND = c1*
			(-0.5 + 1.5/(1.0 + x*x + y*y));
	    }}}
	break;

      default: // this if for patches 0 to 3
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2; 

		    
	      	(*y_a_ptr).u[M_A_T]GRD_IND = c1*
			(-0.5 + 1.5*y*y/(1.0 + x*x + y*y));
		(*y_a_ptr).u[M_B_T]GRD_IND = c1*
			(-0.5 + 1.5*y*y/(1.0 + x*x + y*y));
		(*y_a_ptr).u[M_C_T]GRD_IND = c1*
			(-0.5 + 1.5*y*y/(1.0 + x*x + y*y));
	    }}}
	break; 
      }
    }

      break;
    

    case 2: {// give data using cross-like coordinates
      {register int grid_ind1, grid_ind2;
      FLOAT x=0.;
      FLOAT y=0.;
      FLOAT z=0.;
      FLOAT r = 0.;

      for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
	    x =  x_i+(FLOAT)(grid_ind1)*one_dN1;
	    y =  y_i+(FLOAT)(grid_ind2)*one_dN2; 
	    r = (x-c1_1)*(x-c1_1) + (y-c1_2)*(y-c1_2);

	    if(r >= sigma_b0){
	    }

	    else{/*
	      (*y_a_ptr).u[PHI_T]GRD_IND = c1*
		  pow(r-sigma_b1,8)/pow(sigma_b1,8)    
		+ 
		v1*(*y_a_ptr).u[PHI_X]GRD_IND
		+ 
		v2*(*y_a_ptr).u[PHI_Y]GRD_IND;*/
	    }
	  }}}
    }
      break;
    


/* -------------------------------------------------------------*/
/* -------------------------------------------------------------*/
/* CALSE l=2,m-1 multipole */
/* -------------------------------------------------------------*/
/* -------------------------------------------------------------*/
  case 3: // l=2,m=1 multipole
      switch ((*grid_ptr).grid){ // patches 4 and 5 are the ones that have the poles
      case 0:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
				    /*
	      (*y_a_ptr).u[PHI_T]GRD_IND = c1*(y/(1.0+x*x+y*y));*/
	    }}}
	break;
				
      case 2:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
			
/*	      (*y_a_ptr).u[PHI_T]GRD_IND = c1*(-y/(1.0+x*x+y*y));*/
	    }}}
	break;
		
      case 1:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
				    
/*	      (*y_a_ptr).u[PHI_T]GRD_IND = c1*(-y*x/(1.0+x*x+y*y));*/
	    }}}
	break;
				
      case 3:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
			
/*	      (*y_a_ptr).u[PHI_T]GRD_IND = c1*((y*x)/(1.0+x*x+y*y));*/
	    }}}
	break;
		
      case 4:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
				    
/*	      (*y_a_ptr).u[PHI_T]GRD_IND = c1*(-y/(1.0+x*x+y*y));*/
	    }}}

	break;
      case 5:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
	
/*	      (*y_a_ptr).u[PHI_T]GRD_IND = c1*(-y/(1.0+x*x+y*y));*/
	    }}}
	break; 
      }
    
    break;


/* -------------------------------------------------------------*/
/* -------------------------------------------------------------*/


    case 4: {// l=4 multipole
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
        FLOAT theta, phi;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
	      theta = theta_cub(x,y,(*grid_ptr).grid) ;
              phi = phi_cub(x,y,(*grid_ptr).grid);

	  (*y_a_ptr).u[M_A_T]GRD_IND =  0.0;
	  (*y_a_ptr).u[M_B_T]GRD_IND =  0.0;
	  (*y_a_ptr).u[M_C_T]GRD_IND =  0.0;
	    }}}
      }

      break;
 
    default: // unknown initial data type
      printf("unknown initial data type = %d \n", initial_data_type);
      exit(1);
      break;

    } 


//-------------------------------------------------------

//----------- the inverse metric ------------------------

//-------------------------------------------------------


{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
				    
#include"S2_inv_metric.c"

	      auxfields.u_aux[G_UP_XX]GRD_IND = g_up_xx;
	      auxfields.u_aux[G_UP_YY]GRD_IND = g_up_yy;
	      auxfields.u_aux[G_UP_XY]GRD_IND = g_up_xy;
//printf("(%f,%f), g_up_xx = %f,  g_up_xy = %f, g_up_yy = %f \n",x,y,g_up_xx,g_up_xy,g_up_yy);

	    }}}

  printf("<LI>Inidat finished </br></LI>\n"); 
}



// auxiliary functions



       FLOAT theta_cub(FLOAT x, FLOAT y, int grid){// This function returns the value of the theta angle in spherical coordinates, starting from the
                                                                  // cubito coordinates. We first compute z = cos(theta) and from it theta.
        FLOAT z;
	switch (grid){ // patches 4 and 5 are the ones that have the poles 4 is the north pole

	      	case 0: case 1: case 2: case 3:

              	z = y/sqrt(1.+ x*x + y*y);

		break;		

      		case 4: //top patch

              		z = 1./sqrt(1.+x*x+y*y);
	     
		break;

	        case 5:

	                z = -1./sqrt(1.+x*x+y*y);

		break;

  		default: //unknown grid variable
    		printf("In theta_cub, unknown value for grid %d\n", grid);
    		exit(1);
    		break;
	}
                //                printf("theta_cub = %f \n",theta_cub(x,y,(*grid_ptr).grid));
		//	printf("acos(z) = %f, z = %f \n",acos(z), z);
	return(acos(z));
	}




       FLOAT phi_cub(FLOAT x, FLOAT y, int grid){
	// This function returns the value of the phi angle in spherical coordinates, starting from the
        // cubito coordinates. Near the poles returns +-PI/2
	// tested ok

        FLOAT phi;
	switch (grid){ // patches 4 and 5 are the ones that have the poles 4 is the north pole

	      	case 0: 

              	//phi = asin(2.*x/(1.+x*x))/2.;
                phi = atan(x);
		break;		

		case 2: 

              	phi = atan(x) + PI;

		break;		

		case 1: 

              	phi = atan(x)  + PI/2.;

		break;		

		case 3:

              	phi = atan(x)  + 3.*PI/2.;

		break;		

      		case 4: //top patch

			if(y < -0.00001){
	  	        	phi = -atan(x/y);
			}
			else{
				if(y > 0.00001){
								phi = -atan(x/y) + PI;
								}
								else{if(x>=0){
												phi = PI/2.;
									}
												else{
														phi = -PI/2.;
												}
								}	
				}	
	     
		break;

	        case 5:

			if(y < -0.00001){
	  	        	phi = atan(x/y) + PI;
			}
			else{
				if(y > 0.00001){
	  	        	phi = atan(x/y);
				}
				else{if(x>=0){
	     					phi = PI/2.;
						}
						else{
							phi = -PI/2.;
							}
					}	
			}             		

		break;

  		default: //unknown grid variable
    		printf("In phi_cub, unknown value for grid %d\n", grid);
    		exit(1);
    		break;
	}
	return(phi);
	}









// function to asign different theta profiles
       FLOAT fz(FLOAT z){
			//return(3.*PI*z/2.);
			//return(PI*z*z*z/2.);
                        return(PI*z/2.); 
      		        }

